{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 稀疏矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Scipy` 提供了稀疏矩阵的支持（`scipy.sparse`）。\n",
    "\n",
    "稀疏矩阵主要使用 位置 + 值 的方法来存储矩阵的非零元素，根据存储和使用方式的不同，有如下几种类型的稀疏矩阵：\n",
    "\n",
    "类型|描述\n",
    "---|----\n",
    "`bsr_matrix(arg1[, shape, dtype, copy, blocksize])`\t| Block Sparse Row matrix\n",
    "`coo_matrix(arg1[, shape, dtype, copy])`\t| A sparse matrix in COOrdinate format.\n",
    "`csc_matrix(arg1[, shape, dtype, copy])`\t| Compressed Sparse Column matrix\n",
    "`csr_matrix(arg1[, shape, dtype, copy])`\t| Compressed Sparse Row matrix\n",
    "`dia_matrix(arg1[, shape, dtype, copy])`\t| Sparse matrix with DIAgonal storage\n",
    "`dok_matrix(arg1[, shape, dtype, copy])`\t| Dictionary Of Keys based sparse matrix.\n",
    "`lil_matrix(arg1[, shape, dtype, copy])`\t| Row-based linked list sparse matrix\n",
    "\n",
    "在这些存储格式中：\n",
    "\n",
    "- COO 格式在构建矩阵时比较高效\n",
    "- CSC 和 CSR 格式在乘法计算时比较高效"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 构建稀疏矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "创建一个空的稀疏矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<2x3 sparse matrix of type '<type 'numpy.float64'>'\n",
       "\twith 0 stored elements in COOrdinate format>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coo_matrix((2,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "也可以使用一个已有的矩阵或数组或列表中创建新矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 0)\t1\n",
      "  (0, 1)\t2\n",
      "  (1, 2)\t3\n",
      "  (2, 0)\t4\n",
      "  (2, 2)\t5\n"
     ]
    }
   ],
   "source": [
    "A = coo_matrix([[1,2,0],[0,0,3],[4,0,5]])\n",
    "print A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "不同格式的稀疏矩阵可以相互转化："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "scipy.sparse.coo.coo_matrix"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "scipy.sparse.csr.csr_matrix"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = A.tocsr()\n",
    "type(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以转化为普通矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1, 2, 0],\n",
       "        [0, 0, 3],\n",
       "        [4, 0, 5]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = A.todense()\n",
    "C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "与向量的乘法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1, -3, -1])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = np.array([1,0,-1])\n",
    "A.dot(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还可以传入一个 `(data, (row, col))` 的元组来构建稀疏矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "I = np.array([0,3,1,0])\n",
    "J = np.array([0,3,1,2])\n",
    "V = np.array([4,5,7,9])\n",
    "A = coo_matrix((V,(I,J)),shape=(4,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 0)\t4\n",
      "  (3, 3)\t5\n",
      "  (1, 1)\t7\n",
      "  (0, 2)\t9\n"
     ]
    }
   ],
   "source": [
    "print A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "COO 格式的稀疏矩阵在构建的时候只是简单的将坐标和值加到后面，对于重复的坐标不进行处理："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 0)\t1\n",
      "  (0, 2)\t1\n",
      "  (1, 1)\t1\n",
      "  (3, 3)\t1\n",
      "  (1, 1)\t1\n",
      "  (0, 0)\t1\n",
      "  (0, 0)\t1\n"
     ]
    }
   ],
   "source": [
    "I = np.array([0,0,1,3,1,0,0])\n",
    "J = np.array([0,2,1,3,1,0,0])\n",
    "V = np.array([1,1,1,1,1,1,1])\n",
    "B = coo_matrix((V,(I,J)),shape=(4,4))\n",
    "print B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "转换成 CSR 格式会自动将相同坐标的值合并："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 0)\t3\n",
      "  (0, 2)\t1\n",
      "  (1, 1)\t2\n",
      "  (3, 3)\t1\n"
     ]
    }
   ],
   "source": [
    "C = B.tocsr()\n",
    "print C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 求解微分方程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import lil_matrix\n",
    "from scipy.sparse.linalg import spsolve\n",
    "from numpy.linalg import solve, norm\n",
    "from numpy.random import rand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构建 `1000 x 1000` 的稀疏矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = lil_matrix((1000, 1000))\n",
    "A[0, :100] = rand(100)\n",
    "A[1, 100:200] = A[0, :100]\n",
    "A.setdiag(rand(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "转化为 CSR 之后，用 `spsolve` 求解 $Ax=b$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = A.tocsr()\n",
    "b = rand(1000)\n",
    "x = spsolve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "转化成正常数组之后求解："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_ = solve(A.toarray(), b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "查看误差："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.4310987107687431e-13"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "err = norm(x-x_)\n",
    "err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## sparse.find 函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "返回一个三元组，表示稀疏矩阵中非零元素的 `(row, col, value)`："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 1 3] [0 2 1 3] [3 1 2 1]\n"
     ]
    }
   ],
   "source": [
    "from scipy import sparse\n",
    "\n",
    "row, col, val = sparse.find(C)\n",
    "print row, col, val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## sparse.issparse 函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "查看一个对象是否为稀疏矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse.issparse(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "或者"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse.isspmatrix(B.todense())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还可以查询是否为指定格式的稀疏矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse.isspmatrix_coo(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse.isspmatrix_csr(B)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
